Στατιστική και Οικολογική Μοντελοποίηση

1 Βασικές προϋποθέσεις

Εννοείται ότι πριν ξεκινήσουμε να κάνουμε το οτιδήποτε, έχουμε δημιουργήσει ένα νέο project στο R-Studio, το όνομα του οποίου - για λόγους τους οποίους θα καταλάβετε αργότερα - ΔΕΝ πρέπει να περιέχει:
1. ελληνικούς χαρακτήρες
2. κενά (αντικαταστήστε τα κενά με _)

Το ίδιο ισχύει ΚΑΙ για το file path του εν λόγω φακέλου (π.χ., όχι ‘E:/Ειδικά Μαθήματα Βοτανικής/SDM’, αλλά ’E:/Eidika_Mathimata_Botanikis/SDM).

Ένα άλλο κρίσιμο σημείο είναι το εξής: θα χρειαστεί να εγκαταστήσετε τα πακέτα που πρόκειται να χρησιμοποιήσετε σε αυτό το tutorial1.

2 Εγκατάσταση και φόρτωση βιβλιοθηκών

Όπως κάθε φορά, θα χρειαστεί να εγκαταστήσουμε και να φορτώσουμε τις απαραίτητες βιβλιοθήκες.

## ===========================================================================##
## Install the main packages
## ===========================================================================##
install.packages(c("tidyverse", "rgbif", "raster", "data.table", "rgeos", "rgdal", 
    "gridExtra", "dismo", "biomod2", "sf", "usdm", "speciesgeocodeR", "pacman", 
    "spThin", "CoordinateCleaner", "ggspatial", "scales", "scrubr", "mapr", 
    "rasterVis", "ade4", "viridis", "biogeo"), dependencies = T)
## ===========================================================================##

Ας τις φορτώσουμε:

## ===========================================================================##
## Load them
## ===========================================================================##
pacman::p_load(tidyverse, rgbif, raster, usdm, rgeos, rgdal, usdm, biomod2, 
    gridExtra, countrycode, devtools, sf, spThin, biogeo, speciesgeocodeR, pacman, 
    dismo, CoordinateCleaner, ggspatial, scales, scrubr, mapr, rasterVis, data.table, 
    ade4, viridis)
## ===========================================================================##

3 Διανυσματικά δεδομένα

3.1 Σύνορα χωρών

Πρώτα θα κατεβάσουμε τις γεωγραφικές μεταβλητές για την χώρα που μας ενδιαφέρει (η Ελλάδα στην προκειμένη περίπτωση). Θα πρέπει να εισάγουμε (στην εντολή country στον κώδικα που παρατίθεται πιο κάτω) τον 3-ψήφιο ISO κωδικό για την χώρα μας που τον βρίσκουμε από αυτόν τον ιστότοπο ή με την εντολή ISO_countries <- getData("ISO3") %>% as.data.frame2. Εάν θέλουμε να περιορίσουμε την περιοχή μελέτης σε μια μικρότερη περιοχή, τότε θα πρέπει να χρησιμοποιήσουμε ένα άλλο λογισμικό, το Qgis, ώστε να “κόψουμε” την περιοχή που μας ενδιαφέρει. Προς το παρόν, δεν θα χρειαστεί να το κάνουμε αυτό, αλλά θα περιοριστούμε να κατεβάσουμε τα πλήρη δεδομένα για την Ελλάδα: .

## ===========================================================================##
## Download the data for Greece
## ===========================================================================##
## There are three levels available, you are free to explore them
Greece <- getData("GADM", country = "GRC", level = 3)
## ===========================================================================##

Μπορούμε να αναπαραστήσουμε γραφικά τα δεδομένα αυτά.

## ===========================================================================##
## Plot Greece
## ===========================================================================##
plot(Greece)

Είναι αρκετά εύκολο στην προκειμένη περίπτωση να περιορίσουμε την έκταση του πολυγώνου στην Κρήτη, με τις ακόλουθες εντολές:

## ===========================================================================##
## First we need to see what is inside the object we created
## ===========================================================================##
Greece
## class       : SpatialPolygonsDataFrame 
## features    : 326 
## extent      : 19.37236, 29.6457, 34.80069, 41.74801  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 16
## names       : GID_0, NAME_0,   GID_1,                      NAME_1,                                                                                        NL_NAME_1,     GID_2,         NAME_2,                        NL_NAME_2,       GID_3,   NAME_3,          VARNAME_3,                                                     NL_NAME_3, TYPE_3,    ENGTYPE_3, CC_3, ... 
## min values  :   GRC, Greece, GRC.1_1,                      Aegean,                                                                 <U+0386><U+03B8><U+03C9><U+03C2>, GRC.1.1_1,          Athos, <U+0386><U+03B8><U+03C9><U+03C2>, GRC.1.1.1_1,   Abdera, Abelokipi-Menemeni,                              <U+0386><U+03B8><U+03C9><U+03C2>,  Dímos, Municipality,   NA, ... 
## max values  :   GRC, Greece, GRC.8_1, Thessaly and Central Greece, Tessa<U+03BB><U+03AF>a <U+03BA>a<U+03B9> Ste<U+03C1>e<U+03AC> <U+0395><U+03BB><U+03BB><U+03AC>da, GRC.8.2_1, West Macedonia,           Tessa<U+03BB><U+03AF>a, GRC.8.2.9_1, Zografou,           Zografos, Tessa<U+03BB><U+03BF><U+03BD><U+03AF><U+03BA><U+03B7><U+03C2>,  Dímos, Municipality,   NA, ...
names(Greece)
##  [1] "GID_0"     "NAME_0"    "GID_1"     "NAME_1"    "NL_NAME_1"
##  [6] "GID_2"     "NAME_2"    "NL_NAME_2" "GID_3"     "NAME_3"   
## [11] "VARNAME_3" "NL_NAME_3" "TYPE_3"    "ENGTYPE_3" "CC_3"     
## [16] "HASC_3"
Crete <- subset(Greece, NAME_1 == "Crete")
## ===========================================================================##
## Plot Crete
## ===========================================================================##
plot(Crete)

Κάντε το ίδιο για την Αττική, τον Άθω και την Θεσσαλία.

3.2 GBIF

Στο σημείο αυτό, είμαστε πλέον στην ευχάριστη θέση να μπορούμε να χρησιμοποιήσουμε τα δεδομένα θέσης για τα είδη τα οποία μας ενδιαφέρουν. Υπάρχουν δύο περιπτώσεις, όσον αφορά την απόκτηση τέτοιων δεδομένων:
- να τα έχουμε αποκτήσει πρωτογενώς (δηλαδή να έχουμε εργαστεί στο πεδίο και να έχουμε εντοπίσεις τις θέσεις εμφάνισης του εκάστοτε είδους και των υποπληθυσμών του)
- να τα έχουμε αποκτήσει δευτερογενώς (δηλαδή είτε να έχουμε χρησιμοποιήσει τις διαδικτυακές βάσεις ψηφιοποιημένων - και ενίοτε γεωαναφερμένων - δεδομένων GBIF και BIEN, είτε να έχουμε ψηφιοποιήσει χάρτες εμφάνισης θέσης3)

Το ιδανικό βεβαίως είναι να έχουμε αποκτήσει πρωτογενώς τα στοιχεία αυτά, καθότι στην περίπτωση αυτή, μπορούμε να είμαστε σίγουροι για τις γεωγραφικές συντεταγμένες των υποπληθυσμών του προς μελέτη είδους. Στην πλειονότητα των περιπτώσεων όμως, δεν έχουμε πρόσβαση σε τέτοια δεδομένα, οπότε αναγκαστικά καταφεύγουμε στην δεύτερη λύση και δη στις διαδικτυακές βάσεις δεδομένων. Ευτυχώς για εμάς, υπάρχουν κάποια πακέτα4 στην R, όπως το rgbif και το rbien, τα οποία αυτοματοποιούν την διαδικασία μεταφόρτωσης δεδομένων από τις διαδικτυακές αυτές βάσεις5.

Εμείς θα εργαστούμε με το πακέτο rgbif προς το παρόν.

Σήμερα – όπως και στην προηγούμενη εργαστηριακή άσκηση – θα ασχοληθούμε ένα εκ των έξι μονοτυπικών γενών της ελληνικής χλωρίδας, την Petromarula pinnata.

Fig 1. The endemic species Petromarula pinnata from Crete.

Fig 1. The endemic species Petromarula pinnata from Crete.

Πρώτα χρειάζεται να βεβαιωθούμε ότι το προς μελέτη είδος, υπάρχει όντως εντός της διαδικτυακής βάσης την οποία επιθυμούμε να χρησιμοποιήσουμε, προκειμένου να εξάγουμε δεδομένα. Η εντολή name_suggest() ερευνά όλα τα είδη (taxa καλύτερα) τα οποία μοιάζουν με το είδος το οποίο ψάχνουμε και κρατά μόνο το όνομα των taxa το οποίο ταιριάζει ακριβώς με αυτό το οποίο μας ενδιαφέρει6.

spp_petromarula <- name_suggest(q = "Petromarula pinnata", rank = "species", 
    limit = 10000)
(spp_petromarula <- spp_petromarula[grepl("^Petromarula pinnata$", spp_petromarula$canonicalName), 
    ])
key canonicalName rank
3167989 Petromarula pinnata SPECIES

Στο βήμα αυτό αναγνωρίσαμε τον μοναδικό αριθμό αναφοράς (key) του προς μελέτη είδους και μόνο χρησιμοποιώντας τον αριθμό αυτό μπορούμε να είμαστε βέβαιοι ότι τα δεδομένα θέσης τα οποία θα λάβουμε από την διαδικτυακή αυτή βάση δεδομένων θα αναφέρονται στο είδος το οποίο μας ενδιαφέρει.

Στη συνέχεια, θα δούμε πόσα γεωαναφερμένα δεδομένα θέσης έχουμε στην διάθεση μας.

3.3 Μεταφόρτωση των γεωαναφερμένων δεδομένων θέσης

## ===========================================================================##
## Download data
## ===========================================================================##
sp_occ <- occ_search(taxonKey = spp_petromarula$key, country = "GR", hasCoordinate = T, 
    limit = 1000, return = "data")
## ===========================================================================##


## ===========================================================================##
## To prevent any problem with the pathway it is a good practice to remove
## blank spaces from species names
## ===========================================================================##
sp_occ$name <- sub(" ", "_", sp_occ$name)
## ===========================================================================##


## ===========================================================================##
## Total number of occurrences
## ===========================================================================##
sort(table(sp_occ$name), decreasing = T)
## Petromarula_pinnata A.DC. 
##                        77

Θα πρέπει να τονίσουμε το εξής: Εφ’όσον δεν έχουμε συλλέξει πρωτογενώς τα δεδομένα μας, θα χρειαστεί να ελέγξουμε εάν οι θέσεις εμφάνισης που μεταφορτώσαμε από τις διαδικτυακές βάσεις δεδομένων, συμφωνούν με τα γνωστά όρια εξάπλωσης του είδους που μας ενδιαφέρει. Εάν όντως υπάρχουν κάποια σημεία τα οποία βρίσκονται εκτός της περιοχής εξάπλωσης, θα χρειαστεί να τα αφαιρέσουμε από την ανάλυση.

## ===========================================================================##
## Clean the coordinates
## ===========================================================================##
flags_gbif <- clean_coordinates(x = sp_occ, lon = "decimalLongitude", lat = "decimalLatitude", 
    countries = "countryCode", species = "species", tests = c("capitals", "centroids", 
        "equal", "gbif", "zeros", "seas"), seas_ref = Crete)
## ===========================================================================##


## ===========================================================================##
## Keep only those that are OK
## ===========================================================================##
flags_gbif <- flags_gbif %>% mutate(Data = ifelse(.summary == "TRUE", "Clean", 
    "Flagged"), NAME_1 = "Crete", Dataset = "GBIF", longitude = decimalLongitude, 
    latitude = decimalLatitude, scientific_name = "Petromarula pinnata") %>% 
    dplyr::filter(str_detect(Data, "Cle")) %>% dplyr::select(scientific_name, 
    longitude, latitude, Dataset)
## ===========================================================================##


## ===========================================================================##
## Create a duplicate for safekeeping
## ===========================================================================##
flags_gbif_dpl <- flags_gbif
## ===========================================================================##

Ας φτιάξουμε έναν φάκελο ώστε να αποθηκεύσουμε τα δεδομένα μας και ας τα σώσουμε, προκειμένου να μπορέσουμε να τα φορτώσουμε κάποια άλλη στιγμή:

dir.create("./Species", recursive = TRUE, showWarnings = FALSE)
dir.create("./Species/Petromarula pinnata", recursive = TRUE, showWarnings = FALSE)

saveRDS(flags_gbif, "./Species/ Petromarula pinnata.rds")

4 Μεταφόρτωση περιβαλλοντικών μεταβλητών

Τα κλιματικά δεδομένα μας μπορούμε να τα κατεβάσουμε από το WorldClim, ενώ τα υψομετρικά δεδομένα από εδώ και τα δεδομένα σχετικά με την ξηρασία και την εξατμισιοδιαπνοή από εδώ. - Εάν φυσικά δουλεύουμε για μια συγκεκριμένη εργασία (πτυχιακή, μεταπτυχιακή, κτλ) και όχι στα πλαίσια ενός προπτυχιακού μαθήματος επιλογής. Προς το παρόν, θα χρησιμοποιήσουμε κάποιες εντολές για να κατεβάσουμε τα δεδομένα που θέλουμε μέσω της R.

Όπως μπορείτε να διαπιστώσετε, υπάρχουν πάρα πολλά κλιματικά μοντέλα (επί παραδείγματι, CCSM4, HadGEM-2, κτλ) και τέσσερα (4) διαφορετικά κλιματικά σενάρια που έχουν να κάνουν με τη συγκέντρωση των αερίων του θερμοκηπίου στην ατμόσφαιρα εντός του 21ου αιώνα { ενδεικτική βιβλιογραφία ή άλλη ενδεικτική βιβλιογραφία ή ακόμα μια ενδεικτική βιβλιογραφία ή τελευταία ενδεικτική βιβλιογραφία}7.

Συνεπώς, δεν μπορούμε να χρησιμοποιήσουμε μόνο ένα κλιματικό μοντέλο, πόσω μάλλον ένα κλιματικό σενάριο. Στο σημείο αυτό, είναι εύλογο να αρχίσετε να πανικοβάλεστε8. Ευτυχώς όμως, δεν χρειάζεται και ούτε είναι επιστημονικά σωστό, να χρησιμοποιήσουμε όλα τα κλιματικά μοντέλα. Αναλόγως της περιοχής στην οποία εργαζόμαστε, κάποια κλιματικά μοντέλα αναπαριστούν καλύτερα τις περιβαλλοντικές συνθήκες σε σχέση με κάποια άλλα (McSweeney et al., 2015). Καθώς θα διαπιστώσετε (όταν με το καλό διαβάσετε την προαναφερθείσα εργασία), ακόμα και έτσι, είναι πολλά τα κλιματικά μοντέλα για την περιοχή μας. Όμως, μόνο 3 από αυτά έχουν δεδομένα και για τα 4 κλιματικά σενάρια (BCC, HadGEM, CCSM4). Εμείς χάριν συντομίας, θα εργαστούμε όμως μόνο με ένα (1) κλιματικό μοντέλο και ενα (1) κλιματικό σενάριο.

## ============================================================================##
## First, let's download current climate data for the whole planet
## ============================================================================##
current_climate <- getData("worldclim", var = "bio", res = 2.5)
## ============================================================================##


## ============================================================================##
## Then, let's download future climate data
## ============================================================================##
future_climate <- getData("CMIP5", var = "bio", res = 2.5, rcp = 85, model = "CC", 
    year = 70)

## ============================================================================##

Ας αλλάξουμε τα ονόματα των περιβαλλοντικών μεταβλητών.

names(current_climate) <- c("bio_1", "bio_2", "bio_3", "bio_4", "bio_5", "bio_6", 
    "bio_7", "bio_8", "bio_9", "bio_10", "bio_11", "bio_12", "bio_13", "bio_14", 
    "bio_15", "bio_16", "bio_17", "bio_18", "bio_19")

names(future_climate) <- c("bio_1", "bio_2", "bio_3", "bio_4", "bio_5", "bio_6", 
    "bio_7", "bio_8", "bio_9", "bio_10", "bio_11", "bio_12", "bio_13", "bio_14", 
    "bio_15", "bio_16", "bio_17", "bio_18", "bio_19")

Καθώς τα κλιματικά μας δεδομένα όπως διαπιστώσατε έχουν εξαιρετικά μεγάλο εύρος τιμών, θα χρειαστεί να τα μετατρέψουμε σε oC όσον αφορά την θερμοκρασία (bio_1). Εδώ μπορείτε να διαπιστώσετε το γιατί9. Αυτό γίνεται με την ακόλουθη εντολή:

gain(current_climate$bio_1) = 0.1

gain(future_climate$bio_1) = 0.1

Στη συνέχεια, θα κατεβάσουμε τα υψομετρικά δεδομένα για την χώρα που μας ενδιαφέρει:

altitude <- getData("alt", country = "GRC", mask = TRUE)

4.1 Κοπή των περιβαλλοντικών δεδομένων στην έκταση της Ελλάδας

Στη συνέχεια, θα ‘κόψουμε’ τις περιβαλλοντικές μας μεταβλητές στο σχήμα της Ελλάδας και ακολούθως θα δημιουργήσουμε ένα νέο αρχείο που θα περιέχει τις μεταβλητές αυτές και την περιοχή μελέτης μας:

## ============================================================================##
## First for Greece
## ============================================================================##
current_climate_crop <- crop(current_climate, Greece, snap = "in") %>% mask(., 
    Greece)


future_climate_crop <- crop(future_climate, Greece, snap = "in") %>% mask(., 
    Greece)

## ============================================================================##



## ============================================================================##
## Then for Crete
## ============================================================================##
current_climate_Crete <- crop(current_climate_crop, Crete, snap = "in") %>% 
    mask(., Crete)

future_climate_Crete <- crop(future_climate_crop, Crete, snap = "in") %>% mask(., 
    Crete)

## ============================================================================##

Κάντε το ίδιο για την Αττική, τον Άθω και την Θεσσαλία.

4.2 Γραφική απεικόνιση

Μπορούμε να αναπαραστήσουμε γραφικά τα αποτελέσματα ως εξής:

gplot(current_climate_Crete[[1]]) + geom_raster(aes(fill = value)) + scale_fill_gradientn(colours = c("brown", 
    "red", "yellow", "darkgreen", "green")) + coord_equal() + xlab("Longitude") + 
    ylab("Latitude") + ggtitle("Temperature")

histogram(current_climate_Crete[[1]])

levelplot(current_climate_Crete[[1]])

# plot3D(Current_climate_Crete[[1]])

Κάντε το ίδιο για την Αττική, τον Άθω και την Θεσσαλία.

4.3 Βασικά στατιστικά στοιχεία

Στη συνέχεια, θα υπολογίσουμε κάποια βασικά στατιστικά στοιχεία για κάθε αρχείο το οποίο δημιουργήσαμε:

cellStats(current_climate_crop, "quantile")
##      bio_1 bio_2 bio_3   bio_4 bio_5 bio_6 bio_7 bio_8 bio_9 bio_10 bio_11
## 0%     2.4    55    27 4730.00   180   -95   182   -37    94    106    -53
## 25%   12.2    92    33 5988.75   280   -11   262    60   205    207     37
## 50%   14.5   101    34 6516.00   298    12   288    81   229    230     60
##      bio_12 bio_13 bio_14 bio_15 bio_16 bio_17 bio_18 bio_19
## 0%      385  52.00      0     20  138.0      2      2    110
## 25%     542  74.00      8     32  206.0     34     39    184
## 50%     653 103.00     15     50  269.5     65     69    250
##  [ reached getOption("max.print") -- omitted 2 rows ]
cellStats(current_climate_crop, "mean")
##      bio_1      bio_2      bio_3      bio_4      bio_5      bio_6 
##   14.03543   97.27204   34.09507 6455.04935  295.72570   15.26490 
##      bio_7      bio_8      bio_9     bio_10     bio_11     bio_12 
##  280.46080   80.48347  222.14066  224.08342   60.01553  695.86060 
##     bio_13     bio_14     bio_15     bio_16     bio_17     bio_18 
##  110.14041   15.94489   51.10471  296.70353   61.12325   66.01879 
##     bio_19 
##  269.81864
cellStats(current_climate_crop, "sd")
##      bio_1      bio_2      bio_3      bio_4      bio_5      bio_6 
##   2.854903  14.953890   2.226295 671.676911  25.702217  37.634237 
##      bio_7      bio_8      bio_9     bio_10     bio_11     bio_12 
##  35.202266  31.311853  27.093454  26.693768  33.891841 189.062736 
##     bio_13     bio_14     bio_15     bio_16     bio_17     bio_18 
##  41.171401   9.894798  20.722858 109.157622  32.687959  35.129005 
##     bio_19 
## 102.446221
summary(current_climate_crop)
##           bio_1 bio_2 bio_3    bio_4 bio_5 bio_6 bio_7 bio_8 bio_9 bio_10
## Min.        2.4    55    27  4730.00   180   -95   182   -37    94    106
## 1st Qu.    12.2    92    33  5988.75   280   -11   262    60   205    207
## Median     14.5   101    34  6516.00   298    12   288    81   229    230
##         bio_11 bio_12   bio_13 bio_14 bio_15  bio_16 bio_17 bio_18 bio_19
## Min.       -53    385    52.00      0     20   138.0      2      2    110
## 1st Qu.     37    542    74.00      8     32   206.0     34     39    184
## Median      60    653   103.00     15     50   269.5     65     69    250
##  [ reached getOption("max.print") -- omitted 3 rows ]

Κάντε το ίδιο για την Αττική, τον Άθω και την Θεσσαλία.

Μπορούμε να δημιουργήσουμε υποσύνολα των περιοχών που μας ενδιαφέρουν, ώστε να ανταποκρίνονται σε συγκεκριμένες κλιματικές (ή άλλες συνθήκες):

subset_GR <- current_climate_crop[[1]] < 15 & current_climate_crop[[1]] > 10.1

plot(subset_GR)

Δοκιμάστε το και για άλλες περιοχές και μεταβλητές.

4.4 Σύγκριση

Τώρα, θα πρέπει να ελέγξουμε εάν τα υψομετρικά δεδομένα έχουν την ίδια ανάλυση (resolution), μέγεθος (nrow, ncol) και έκταση (bbox), με τα περιβαλλοντικά μας δεδομένα:

altitude
## class      : RasterLayer 
## dimensions : 852, 1260, 1073520  (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333  (x, y)
## extent     : 19.2, 29.7, 34.7, 41.8  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## source     : G:/Academic/Lessons/EMB/Labs/Labs/Fourth_Lab/Fourth_Lab/GRC_msk_alt.grd 
## names      : GRC_msk_alt 
## values     : -11, 2724  (min, max)
current_climate_crop
## class      : RasterBrick 
## dimensions : 165, 246, 40590, 19  (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667  (x, y)
## extent     : 19.375, 29.625, 34.83333, 41.70833  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      :  bio_1,  bio_2,  bio_3,  bio_4,  bio_5,  bio_6,  bio_7,  bio_8,  bio_9, bio_10, bio_11, bio_12, bio_13, bio_14, bio_15, ... 
## min values :    2.4,   55.0,   27.0, 4730.0,  180.0,  -95.0,  182.0,  -37.0,   94.0,  106.0,  -53.0,  385.0,   52.0,    0.0,   20.0, ... 
## max values :   19.3,  123.0,   41.0, 7606.0,  348.0,  104.0,  341.0,  155.0,  270.0,  270.0,  133.0, 1336.0,  226.0,   49.0,  102.0, ...

4.5 Resolution change

Στην περίπτωση κατά την οποία τα δύο σετ δεδομένων έχουν διαφορετική ανάλυση, έκταση και μέγεθος, θα χρειαστεί να αλλάξουμε την ανάλυση των υψομετρικών δεδομένων10.

agg_altitude <- aggregate(altitude, 
                 5, ## Why choose 5 and not any other number?
                 fun = mean)

4.6 Pixel size

Μετά, θα αλλάξουμε το μέγεθος των pixel των υψομετρικών δεδομένων για να συμφωνεί με αυτό των περιβαλλοντικών δεδομένων:

altitude_rsmp <- resample(agg_altitude, current_climate_crop, method = "bilinear")

4.7 Final stack

Τέλος, θα τα ενώσουμε σε ένα αρχείο, τόσο για το παρόν, όσο και για το μέλλον:

## ============================================================================##
## Merge (stack) climate and altitude
## ============================================================================##
current_climate_crop <- stack(current_climate_crop, altitude_rsmp)
future_climate_crop <- stack(future_climate_crop, altitude_rsmp)
## ============================================================================##

4.8 Αποθήκευση σε συμπιεσμένη μορφή

Καλό είναι να αποθηκεύσετε τα αρχεία αυτά υπό την μορφή Rds και να τα αφαιρέσετε από την μνήμη της R11. Αυτό το κάνουμε, προκειμένου να μπορούμε εύκολα να το φορτώνουμε στην R και να μην υπάρχουν τα δεδομένα αυτά συνεχώς στην μνήμη της session που έχουμε ανοικτή, διότι αυξάνουν αρκετά τον χρόνο απόκρισης της (καθυστερεί πολύ η R). Μπορείτε ανά πάσα στιγμή να τα φορτώσετε στην R, με την ακόλουθη εντολή:

## Creating a folder to save our data
dir.create("./Data/", recursive = TRUE, showWarnings = FALSE)

## Save our data in rds format
saveRDS(current_climate_crop, "./Data/Environmental variables and Elevation.rds")
saveRDS(future_climate_crop, "./Data/Environmental variables and Elevation 2070.rds")

## Load our data from rds format
var_GR <- readRDS("./Data/Environmental variables and Elevation.rds")
var_GR_70 <- readRDS("./Data/Environmental variables and Elevation 2070.rds")

4.9 Έλεγχος συγγραμικότητας και Variance Inflation Factor (VIF)

Όπως θα γνωρίζετε, ΔΕΝ μπορούμε να χρησιμοποιήσουμε σε μια ανάλυση δύο μεταβλητές, οι οποίες εμφανίζουν συντελεστή συσχέτισης \(> 0.7\)12, καθότι στην προκειμένη περίπτωση οι δύο αυτές μεταβλητές είναι συγγραμικές και ως εκ τούτου εάν τις συμπεριλάβουμε και τις δύο στην ανάλυση μας, τότε δεν θα γνωρίζουμε σε ποιά εκ των δύο θα οφείλεται το αποτέλεσμα της ανάλυσης μας.

Προκειμένου να ξεπεράσουμε αυτόν τον σκόπελο, θα τρέξουμε την ακόλουθη εντολή13:
.

uncorrelated_vars <- vifcor(current_climate_Crete %>% as.data.frame(), th = 0.7)
uncorrelated_vars <- uncorrelated_vars@results$Variables[uncorrelated_vars@results$VIF <= 
    5]
uncorrelated_vars <- droplevels(uncorrelated_vars)

current_climate_Crete <- stack(subset(current_climate_Crete, levels(uncorrelated_vars)))

future_climate_Crete <- stack(subset(future_climate_Crete, levels(uncorrelated_vars)))

5 Δεδομένα θέσης και περιβαλλοντικές μεταβλητές

Πλέον μπορούμε να συνδέσουμε τις γεωαναφερμένες θέσεις εμφάνισης του είδους που μας ενδιαφέρει με δεδομένα ψηφιδοπλέγματος. Στο σημείο αυτό, μπορούμε να εξάγουμε τα περιβαλλοντικά δεδομένα για κάθε μια εκ των θέσεων εμφάνισης του είδους που μας ενδιαφέρει με μια μόνο εντολή, αναδεικνύοντας κατ’αυτόν τον τρόπο το πλεονέκτημα του να έχουμε ενώσει προηγουμένως τα raster αρχεία μας σε ένα RasterStack.

coordinates(flags_gbif) <- c("longitude", "latitude")

pts.clim <- raster::extract(current_climate_Crete, flags_gbif, method = "bilinear")

sp_occ_clim <- data.frame(cbind(coordinates(flags_gbif), pts.clim, flags_gbif@data))

Και τώρα, ήρθε η στιγμή να οπτικοποιήσουμε τα αποτελέσματα μας και δη, την παρουσία του είδους μας στην περιοχή μελέτης:

map.ext <- extent(Crete)

plot(current_climate_crop[[20]], col = terrain.colors(100), ext = map.ext)

plot(flags_gbif, pch = 16, cex = 0.5, add = T)

plot(Crete, add = T)

Στο σημείο αυτό, θα ήθελα να σας τονίσω ότι ορισμένες φορές, πέρα από τους ελέγχους που διενεργήσαμε έως τώρα, ίσως χρειαστεί να προχωρήσουμε σε μια αραίωση (thinning) των θέσεων εμφάνισης. Αυτό μπορεί να είναι απαραίτητο όταν έχουμε πάρα πολλά σημεία τα οποία είναι αρκετά κοντά το ένα στο άλλο14. Σε κάθε περίπτωση, το επιστημονικά ορθό είναι να υπάρχει μια (1) θέση εμφάνισης ανά km2. Ένα πακέτο που μπορεί να μας βοηθήσει είναι το spThin15.

## Min. distance: 1000 m Number of reps: 1000
spoccs <- flags_gbif_dpl %>% as.data.frame() %>% dplyr::rename(x = longitude, 
    y = latitude, Species = scientific_name) %>% addmainfields(species = "Species")

thinned <- thin(loc.data = spoccs, lat.col = "y", long.col = "x", spec.col = "Species", 
    thin.par = 5, reps = 1000, verbose = T, out.dir = "./Thinned", locs.thinned.list.return = TRUE, 
    write.files = TRUE, out.base = "Petromarula pinnata thinned")
## ********************************************** 
##  Beginning Spatial Thinning.
##  Script Started at: Sat May 23 12:17:05 2020
## lat.long.thin.count
##   39 
## 1000 
## [1] "Maximum number of records after thinning: 39"
## [1] "Number of data.frames with max records: 1000"
## [1] "Writing new *.csv files"
## [1] "Writing file: ./Thinned/Petromarula pinnata thinned_thin1_new.csv"
## [1] "Writing file: ./Thinned/Petromarula pinnata thinned_thin2_new.csv"
## [1] "Writing file: ./Thinned/Petromarula pinnata thinned_thin3_new.csv"
## [1] "Writing file: ./Thinned/Petromarula pinnata thinned_thin4_new.csv"
## [1] "Writing file: ./Thinned/Petromarula pinnata thinned_thin5_new.csv"

6 Μοντελοποίηση

Όλα τα απαιτούμενα δεδομένα και αρχεία θα τα βρείτε στον ιστότοπο του μαθήματος στο εδώ

6.1 Εισαγωγή

Είμαστε πλεόν στην ευχάριστη θέση να έχουμε στην διάθεση μας όλα τα δεδομένα τα οποία είναι απαραίτητα για να προχωρήσουμε στο κυρίως σκέλος της ανάλυσης μας. Κανονικά θα έπρεπε να χρησιμοποιήσουμε πολλά μοντέλα, τα οποία ανήκουν στην ίδια ομάδα (Παρουσίας/Απουσίας ή Ψεύδο-απουσίας - Presence/Absence: P/A), αλλά και σε άλλη ομάδα (Μόνο Παρουσίας - Presence Only: P/O). Κάθε ένα εκ των μοντέλων αυτών έχει τα πλεονεκτήματα και τα μειονεκτήματα του16. Στο σημείο αυτό, πρέπει να τονίσουμε ότι είναι εντελώς λάθος να χρησιμοποιούμε δύο μοντέλα17 τα οποία ανήκουν σε διαφορετικές ομάδες (εν αντιθέσει με ότι προτείνεται στο παράδειγμα του προαναφερθέντος βιβλίου), αλλά να χρησιμοποιεί πολλά μοντέλα τα οποία ανήκουν στην ίδια ομάδα (π.χ., P/A), προκειμένου να μπορεί να βγάλει ασφαλή συμπεράσματα σχετικά με το ερώτημα το οποίο καλείται να απάντησει και να είναι σίγουρος ότι δεν πρόκειται περί στατιστικών κατασκευασμάτων (Araújo & New, 2007).

6.2 Φόρτωση των δεδομένων θέσης

Τώρα θα φορτώσουμε ένα από τα set δεδομένων18 τα οποία έχουν προέλθει από το thinning των θέσεων εμφάνισης19.

thinned_1 <- readxl::read_excel("./Thinned/Petromarula pinnata thinned_thin1.xlsx")

6.3 Μορφοποιήση των δεδομένων

Το πρώτο πράγμα το οποίο θα χρειαστεί να κάνουμε, είναι να μορφοποιήσουμε με τον κατάλληλο τρόπο τα δεδομένα μας, ώστε να μπορέσουμε να τρέξουμε τις εντολές από τα πακέτα biomod και ecospat. Το δεύτερο πακέτο θα το χρησιμοποιήσουμε στην περίπτωση όπου οι θέσεις εμφάνισης του είδους που μας ενδιαφέρει είναι \(<50\)20.

Ένα άλλο σημείο το οποιο πρέπει να προσέξουμε, έχει να κάνει με τον αριθμό των απουσιών (absences - στην πλειονότητα των περιπτώσεων βέβαια πρόκειται για ψευδοαπουσίες -pseudoabsences, καθώς στην πραγματικότητα δεν έχουμε δεδομένα απουσίας του εκάστοτε είδους). Ο αριθμός αυτός πρέπει να είναι ίσος με τον αριθμό των παρουσιών21.

Τέλος, ένα τρίτο, κρίσιμο σημείο είναι ο αριθμός των set των pseudoabsences που θα περιλαμβάνει το μοντέλο μας, προκειμένου να αποφύγουμε τον σκόπελο της μεροληπτικής δειγματοληψίας (sampling bias)22.

sp_formated_data <- BIOMOD_FormatingData(resp.var = rep(1, nrow( thinned_1 ) ),
                                         
                                    ## Our rasterStack with the uncorrelated
                                    ## variables
                                    expl.var = current_climate_Crete,
                                    
                                    ## The names of our coords
                                    resp.xy = thinned_1[,c('x', 'y')],
                                    
                                    ## Species' name
                                    resp.name = "Petromarula pinnata",
                                    
                                    ## Number of pseudoabsences' sets
                                    PA.nb.rep = 2,
                                    
                                    ## Number equaling that of the presences
                                    PA.nb.absences = nrow(thinned_1), 
                                    
                                    ## You can change this to 'disk' if you'd
                                    ## like
                                    PA.strategy = 'random',
                                    
                                    ## Min distance from occurrences
                                    PA.dist.min  =  NULL,
                                    
                                    ## Max distance from occurrences
                                    PA.dist.max  =  NULL) 
## 
## -=-=-=-=-=-=-=-=-=-= Petromarula pinnata Data Formating -=-=-=-=-=-=-=-=-=-=
## 
##  Response variable name was converted into Petromarula.pinnata
##       ! No data has been set aside for modeling evaluation
##    > Pseudo Absences Selection checkings...
##    > random pseudo absences selection
##    > Pseudo absences are selected in explanatory variables
##          ! Some NAs have been automaticly removed from your data
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= Done -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
sp_formated_data
## 
## -=-=-=-=-=-=-=-=-=-=-=-= 'BIOMOD.formated.data.PA' -=-=-=-=-=-=-=-=-=-=-=-=
## 
## sp.name =  Petromarula.pinnata
## 
##   31 presences,  0 true absences and  76 undifined points in dataset
## 
## 
##   4 explanatory variables
## 
##      bio_13          bio_15          bio_3           bio_5      
##  Min.   :106.0   Min.   :77.00   Min.   :29.00   Min.   :208.0  
##  1st Qu.:136.0   1st Qu.:83.00   1st Qu.:30.00   1st Qu.:270.0  
##  Median :155.0   Median :84.00   Median :31.00   Median :277.0  
##  Mean   :159.5   Mean   :83.85   Mean   :31.26   Mean   :274.2  
##  3rd Qu.:182.5   3rd Qu.:85.00   3rd Qu.:32.00   3rd Qu.:284.0  
##  Max.   :215.0   Max.   :89.00   Max.   :33.00   Max.   :301.0  
## 
## 
##  2 Pseudo Absences dataset available ( PA1 PA2 ) with  39 
## absences in each (true abs + pseudo abs)
## 
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

Ας αναπαραστήσουμε γραφικά ότι φτιάξαμε στο προηγούμενο βήμα:

plot(sp_formated_data)

Η εικόνα αυτή μας δείχνει που βρίσκονται οι ψευδοαπουσίες (pseudoabsences) στο set των ψευδοαπουσιών σε σχέση με τις θέσεις εμφάνισης. Όπως αναμενόταν, οι ψευδοαπουσίες είναι τυχαία κατανεμημένες.

6.4 Παραμετροποιήση των μοντέλων

Πλέον μπορούμε να παραμετροποιήσουμε τα μοντέλα που θέλουμε να τρέξουμε:

modelling_options <- BIOMOD_ModelingOptions()

Τώρα, προσδιορίζουμε ποιά μοντέλα θα τρέξουμε, καθώς και πολλές άλλες λεπτομέρειες:

## Always change the name of this object, so as to match the name of your species
Petromarula_pinnata_models <- BIOMOD_Modeling( data = sp_formated_data,
                                       
                                       
                                       ## Here, we indicate which models we are
                                       ## going to run (remember that they need
                                       ## to be from the same group P/A ή P/O).
                                       ## Here, we are going to use some models
                                       ## belonging to the P/A group
                                       ## ('GLM','GBM','GAM','CTA','ANN','SRE',
                                       ## 'FDA','MARS','RF')
                                       models  =  c('RF', 'CTA'),
                                       
                                       models.options = modelling_options,
                                       
                                       ## How many times we will run the
                                       ## analysis in order to evaluate its'
                                       ## results [min: 10 - usually: 100]
                                       NbRunEval = 1,
                                       
                                       ## Training-Testing data split
                                       DataSplit = 80, 
                                       
                                       ## How many times we will run the
                                       ## analysis in order to evaluate the
                                       ## importance of each IV [min: 10]
                                       VarImport = 1,
                                       
                                       ## Evaluation criteria
                                       models.eval.meth  =  c('TSS','ROC','KAPPA'), 
                                       
                                       ## Whether or not we will save the results
                                       SaveObj  =  TRUE,
                                       
                                       ## Should the models be re-scaled in 0-1000 range?
                                       rescal.all.models  =  TRUE,
                                       
                                       ## Should our models be weighted based on all the data?
                                       do.full.models  =  TRUE, 
                                       
                                       ## The name of our analysis
                                       modeling.id  =  "allmodels")

Ας φτιάξουμε τους φακέλους όπου θα αποθηκεύσουμε τα μοντέλα μας:

## Create some folders to store our results
dir.create("./BIOMOD2 Modelling Outpout/", recursive = TRUE, showWarnings = FALSE)

dir.create("./BIOMOD2 Modelling Outpout/Species/", recursive = TRUE, showWarnings = FALSE)

dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata", recursive = TRUE, 
    showWarnings = FALSE)

dir.create("./BIOMOD2 Modelling Outpout/Species/ Petromarula_pinnata/Current", 
    recursive = TRUE, showWarnings = FALSE)


## Save the results in rds format
saveRDS(Petromarula_pinnata_models, "./BIOMOD2 Modelling Outpout/Species/ Petromarula_pinnata/Current/ Petromarula_pinnata models.rds")

6.5 Αξιολόγηση των αρχικών μοντέλων

Ας δημιουργήσουμε ένα αντικείμενο το οποίο θα περιέχει την αξιολόγηση των μοντέλων:

Petromarula_pinnata_scores <- get_evaluations(Petromarula_pinnata_models)

6.6 Σημαντικότητα μεταβλητών

Ας δούμε τη σημαντικότητα των μεταβλητών για κάθε μοντέλο (RF, CTA), κάθε σετ δεδομένων (PA1-2) και κάθε επανάληψη (Run, Full):

(Petromarula_pinnata_var_import <- get_variables_importance(Petromarula_pinnata_models))

Ας δούμε τη σημαντικότητα των μεταβλητών για κάθε μοντέλο (RF, CTA) συνολικά (δηλαδή τον μέσο όρο):

options(digits = 2)
apply(Petromarula_pinnata_var_import, c(1, 2), mean)
##           RF  CTA
## bio_13 0.529 0.75
## bio_15 0.069 0.00
## bio_3  0.053 0.00
## bio_5  0.371 0.48

6.7 Ensemble modelling

Τώρα θα τρέξουμε την ανάλυση που θα ενώνει τα αποτελέσματα των δύο μοντέλων (ονομάζεται ensemble modelling αυτή η τεχνική ανάλυσης):

Petromarula_pinnata_ensemble_models <- BIOMOD_EnsembleModeling( modeling.output = Petromarula_pinnata_models,
                                                   em.by = 'all',
                                                   chosen.models  =  'all',
                                                   
                                                   # You can change that
                                                   eval.metric = 'TSS', 
                                                   
                                                   # Keep only the models with
                                                   # TSS >= 0.8
                                                   eval.metric.quality.threshold = 0.8, 
                                                   models.eval.meth = c('KAPPA','TSS','ROC'),
                                                   prob.mean  =  T, 
                                                   committee.averaging = TRUE,
                                                   prob.mean.weight  =  T, 
                                                   prob.mean.weight.decay  =  'proportional',
                                                   VarImport = 0 )
dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Ensemble", 
    recursive = TRUE, showWarnings = FALSE)

saveRDS(Petromarula_pinnata_ensemble_models, "./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Ensemble/Petromarula_pinnata ensemble.rds")

6.8 Αξιολόγηση των συνολικών μοντέλων (ensemble modelling)

Ας δημιουργήσουμε ένα αντικείμενο το οποίο θα περιέχει την αξιολόγηση των μοντέλων από το ensemble modelling. Όσο πιο κοντά στη μονάδα είναι οι τιμές για τους δείκτες TSS, ROC και KAPPA, τόσο καλύτερα23:

options(digits = 4)
(Petromarula_pinnata_ensemble_models_scores <- get_evaluations(Petromarula_pinnata_ensemble_models))
## $Petromarula.pinnata_EMmeanByTSS_mergedAlgo_mergedRun_mergedData
##       Testing.data Cutoff Sensitivity Specificity
## KAPPA        0.892  621.0         100       93.42
## TSS          0.934  621.0         100       93.42
## ROC          0.983  618.5         100       93.42
## 
## $Petromarula.pinnata_EMcaByTSS_mergedAlgo_mergedRun_mergedData
##       Testing.data Cutoff Sensitivity Specificity
## KAPPA        0.869    747       96.77       93.42
## TSS          0.902    747       96.77       93.42
## ROC          0.960    750       96.77       93.42
## 
## $Petromarula.pinnata_EMwmeanByTSS_mergedAlgo_mergedRun_mergedData
##       Testing.data Cutoff Sensitivity Specificity
## KAPPA        0.892  621.0         100       93.42
## TSS          0.934  621.0         100       93.42
## ROC          0.983  619.5         100       93.42
## ===========================================================##
## Vizualisation
## ===========================================================##

g2 <- models_scores_graph(Petromarula_pinnata_ensemble_models, by = "models", 
    metrics = c("TSS", "KAPPA"))

6.9 Προβολή στο μέλλον

Ας τρέξουμε την ανάλυση για το πώς θα αλλάξει η κατάσταση του είδους που μας ενδιαφέρει στο μέλλον. Αρχικά θα τρέξουμε την ανάλυση για κάθε μοντέλο ξεχωριστά (BIOMOD_Projection) και συνολικά (BIOMOD_EnsembleForecasting) για τις παρούσες συνθήκες και στη συνέχεια για τις μελλοντικές συνθήκες (με τις ίδιες εντολές, αλλά θα αλλάξουμε το raster αρχείο στην εντολή new.env):

Petromarula_pinnata_models_proj_current <- BIOMOD_Projection(modeling.output = Petromarula_pinnata_models, 
    new.env = current_climate_Crete, proj.name = "current", binary.meth = "TSS", 
    selected.models = "all", build.clamping.mask = TRUE, output.format = ".img", 
    do.stack = T)
Petromarula_pinnata_ensemble_models_proj_current <- BIOMOD_EnsembleForecasting(EM.output = Petromarula_pinnata_ensemble_models, 
    projection.output = Petromarula_pinnata_models_proj_current, binary.meth = "TSS", 
    output.format = ".img", do.stack = T)
dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/", 
    recursive = TRUE, showWarnings = FALSE)

dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/Individual", 
    recursive = TRUE, showWarnings = FALSE)

dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/Ensemble", 
    recursive = TRUE, showWarnings = FALSE)

saveRDS(Petromarula_pinnata_models_proj_current, "./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/Individual/Petromarula_pinnata Current Individual Projections.rds")

saveRDS(Petromarula_pinnata_ensemble_models_proj_current, "./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/Current/Projections/Ensemble/Petromarula_pinnata Current Ensemble Projections.rds")
Petromarula_pinnata_proj_2070_cc85 <- BIOMOD_Projection(modeling.output = Petromarula_pinnata_models, 
    new.env = future_climate_Crete, proj.name = "2070_CC85", binary.meth = "TSS", 
    output.format = ".img", do.stack = T)
dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070", recursive = TRUE, 
    showWarnings = FALSE)

dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/", 
    recursive = TRUE, showWarnings = FALSE)

dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/Individual", 
    recursive = TRUE, showWarnings = FALSE)

dir.create("./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/Ensemble", 
    recursive = TRUE, showWarnings = FALSE)

saveRDS(Petromarula_pinnata_proj_2070_cc85, "./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/Individual/Petromarula_pinnata 2070 Individual Projections.rds")
Petromarula_pinnata_ensemble_models_proj_2070_cc85 <- BIOMOD_EnsembleForecasting(EM.output = Petromarula_pinnata_ensemble_models, 
    projection.output = Petromarula_pinnata_proj_2070_cc85, binary.meth = "TSS", 
    output.format = ".img", do.stack = T)
saveRDS(Petromarula_pinnata_ensemble_models_proj_2070_cc85, "./BIOMOD2 Modelling Outpout/Species/Petromarula_pinnata/2070/Projections/Ensemble/Petromarula_pinnata 2070 Ensemble Projections.rds")

6.10 Οπτικοποιήση των συνολικών μοντέλων

Τέλος, ας οπτικοποιήσουμε το αποτέλεσμα της ανάλυσης μας. Πρώτα θα θα δημιουργήσουμε ένα αρχείο το οποίο θα περιέχει την αξιολόγηση της κατάστασης της Petromarula pinnata μέσω της τεχνικής ensemble modelling. Στη συνέχεια, θα κρατήσουμε μόνο τα raster layers για το mean και το weighted mean. Τέλος, θα δημιουργήσουμε ένα αρχείο το οποίο θα περιέχει μόνο τις συντεταγμένες των θέσεων εμφάνισης του είδους που μας ενδιαφέρει:

stk_Petromarula_pinnata_ef_2070_cc85 <- get_predictions(Petromarula_pinnata_ensemble_models_proj_2070_cc85)

stk_Petromarula_pinnata_ef_2070_cc85 <- subset(stk_Petromarula_pinnata_ef_2070_cc85, 
    grep("EMmean|EMwmean", names(stk_Petromarula_pinnata_ef_2070_cc85)))

names(stk_Petromarula_pinnata_ef_2070_cc85) <- sapply(strsplit(names(stk_Petromarula_pinnata_ef_2070_cc85), 
    "_"), getElement, 2)

Petromarula_pinnata_coords <- thinned_1 %>% dplyr::select(x, y)

Στη συνέχεια θα οπτικοποιήσουμε το αποτέλεσμα της ανάλυσης μας έως αυτό το σημείο:

library(viridis)
levelplot(stk_Petromarula_pinnata_ef_2070_cc85, main = "Petromarula pinnata ensemble projections\nin 2070 with cc85", 
    col.regions = viridis(256)) + layer(sp.points(SpatialPoints(Petromarula_pinnata_coords), 
    pch = 20, col = 1))

6.11 Μεταβολή περιοχής εξάπλωσης

Έχουμε φτάσει πλέον στο τελικό σημείο της ανάλυσης μας, όπου θα φορτώσουμε τα αποτελέσματα της ανάλυσης για το mean και το weighted mean της Petromarula pinnata και στη συνέχεια θα υπολογίσουμε και θα οπτικοποιήσουμε τη μεταβολή της περιοχής εξάπλωσης της στο μέλλον:

##==================================
## Load the the data
##==================================

## Current
Petromarula_pinnata_bin_proj_current <- stack('Petromarula.pinnata/proj_current/proj_current_Petromarula.pinnata_ensemble_TSSbin.img')
names(Petromarula_pinnata_bin_proj_current) <- c('ca', 'm', 'wm')

## 2070
Petromarula_pinnata_bin_proj_2070_CC85 <- stack('Petromarula.pinnata/proj_2070_CC85/proj_2070_CC85_Petromarula.pinnata_ensemble_TSSbin.img')
names(Petromarula_pinnata_bin_proj_2070_CC85) <- c('ca', 'm', 'wm')

##=================================
## Species Range Change
##=================================

SRC_current_2070_CC85 <- BIOMOD_RangeSize( Petromarula_pinnata_bin_proj_current,
                                           Petromarula_pinnata_bin_proj_2070_CC85 )
SRC_current_2070_CC85$Compt.By.Models
##    Loss Stable0 Stable1 Gain PercLoss PercGain SpeciesRangeChange
## ca   96     172      50  148    65.75    101.4              35.62
## m    95     168      52  151    64.63    102.7              38.09
## wm   96     172      51  147    65.31    100.0              34.69
##    CurrentRangeSize FutureRangeSize.NoDisp FutureRangeSize.FullDisp
## ca              146                     50                      198
## m               147                     52                      203
## wm              147                     51                      198
Petromarula_pinnata_src_map <- SRC_current_2070_CC85$Diff.By.Pixel
names(Petromarula_pinnata_src_map) <- c("ca cur-2070", "m cur-2070", "wm cur-2070")

##=================================
## Vizualisation
##=================================

my.at <- seq(-2.5,1.5,1)

myColorkey <- list(
                   ## color change
                   at = my.at, 
                   labels = list(
                     
                     ## labels
                     labels = c("lost", "pres", "abs","gain"),
                     
                     ## legend
                     at = my.at[-1]-0.5 
                   ))

rasterVis::levelplot( SRC_current_2070_CC85$Diff.By.Pixel,
                      main = "Petromarula pinnata range change",
                      colorkey = myColorkey,
                      layout = c(1,3),
                      )

7 Εργασία για το σπίτι

Έχετε διορία δεκατεσσάρων (14) ημερών να στείλετε την εργασία σας σε μορφή PDF24 σε αυτό το e-mail.
1. Κατεβάστε γεωγραφικά δεδομένα για την Ελλάδα.
2. Δημιουργήστε ξεχωριστά αρχεία για την Κρήτη, την Αττική, τον Άθω, την Μακεδονία και τη Θεσσαλία.
3. Κατεβάστε τα απαραίτητα κλιματικά δεδομένα.
4. Κατεβάστε υψομετρικά δεδομένα για την Ελλάδα.
5. Υπολογίστε τις βασικές στατιστικές παραμέτρους για όλα τα αρχεία τα οποία δημιουργήσατε.
6. Υπολογίστε τον συντελεστή συσχέτισης για όλα τα περιβαλλοντικά δεδομένα για όλα τα αρχεία που έχετε δημιουργήσει. Εάν υπάρχουν κάποια ζεύγη ανεξάρτητων μεταβλητών με σ.σ. \(>0.7\), θα τις χρησιμοποιήσετε όλες;
7. Κατεβάστε δεδομένα θέσης για τα είδη Asperula pubescens, Centaurea idaea και Scutellaria sieberi.
8. Παραθέστε πληροφορίες σχετικά με τα είδη αυτά (Καθεστώς κινδύνου/προστασίας, καθεστώς ενδημισμού, ταξινομικές πληροφορίες). Προχωρήστε στις απαραίτητες ενέργειες, έτσι ώστε να διαπιστώσετε εάν τα δεδομένα αυτά εντοπίζονται στη σωστή γεωγραφικά περιοχή.
9. Πόσες θέσεις εμφάνισης αριθμούσαν τα είδη αυτά πριν και μετά την διαδικασία αραίωσης;
10. Φορτώστε τα τελικά δεδομένα θέσης για κάθε είδος για το οποίο έχετε κατεβάσει δεδομένα θέσης.
11. Αλλάξτε τα rasterStack αρχεία με αυτά για την Κρήτη σε όλα τα απαραίτητα σημεία του κώδικα. Αλλάξτε επίσης το όνομα του είδους σας. Αλλάξτε τον αριθμό των σετ σε 10 και τον αριθμό των ψευδοαπουσιών25. Αλλάξτε τον αριθμο των φορών που θα τρέξει η ανάλυση σε 10.
12. Ποιά είναι η σημαντικότερη μεταβλητή σε κάθε αλγόριθμο; Συμπίπτουν μεταξύ των αλγορίθμων;
13. Ποιό είναι το καλύτερο μοντέλο σύμφωνα με τους δείκτες αξιολόγησης TSS και KAPPA;.
14. Προβάλετε τα μοντέλα σας στο μέλλον. Παραθέστε τους σχετικούς χάρτες.
15. Πώς μεταβάλλεται η περιοχή εξάπλωσης των ειδών αυτών στο μέλλον;26 Τι συμβαίνει στην περίπτωση της μηδενικής και της πλήρους διασποράς, αντίστοιχα;
16. Παραθέστε το σχετικό γράφημα.

Βιβλιογραφία

Araújo, M.B. & New, M. (2007) Ensemble forecasting of species distributions. Trends Ecol. Evol., 22, 42–47.

Guisan, A., Thuiller, W., & Zimmermann, N.E. (2017) Habitat Suitability and Distribution Models with Applications in R. Cambridge University Press,

McSweeney, C.F., Jones, R.G., Lee, R.W., & Rowell, D.P. (2015) Selecting CMIP5 GCMs for downscaling over multiple regions. Clim. Dyn., 44, 3237–3260.

  1. Υπάρχουν δύο τρόποι για να το κάνετε αυτό γρήγορα:
    1. library(easypackages) packages('package_name_1', 'package_name_2')
    2. install.packages(c('package_name_1', 'package_name_2'), dependecies = TRUE)↩

  2. Ή με αυτήν την εντολή εάν θέλουμε μόνο την Ελλάδα: ISO_countries <- getData("ISO3") %>% as.data.frame %>% filter(NAME=="Greece")↩

  3. Ένα τυπικό παράδειγμα είναι ο Άτλαντας της Χλωρίδας του Αιγαίου ή ο Άτλαντας της Χλωρίδας της Ευρώπης, όμως στην περίπτωση αυτή υπεισέρχονται άλλα προβλήματα τα οποία σχετίζονται με την ανάλυση των εν λόγω χαρτών (π.χ., ο Άτλαντας της Χλωρίδας της Ευρώπης είναι σε ανάλυση \(50 \times 50 km^2\))↩

  4. Δεν είναι τα μοναδικά όμως: υπάρχουν πολλά άλλα, όπως το dismo και το spocc↩

  5. Κάθε ένα από αυτά τα πακέτα, έχει τουλάχιστον μια σύντομη περιγραφή (αγγλιστί: vignette), την οποία είναι ΑΠΑΡΑΙΤΗΤΟ να διαβάσετε - ειδικά αυτές του rgbif↩

  6. Με την εντολή ?rgbif::name_suggest() μπορείτε να δείτε τα arguments που δέχεται η εντολή αυτή και τι μπορείτε να αλλάξετε - σας συνιστώ εντόνως να το κάνετε, αλλά και για κάθε εντολή που χρησιμοποιείτε↩

  7. Καλό είναι να τις διαβάσετε, γιατί θα σας βοηθήσει και στην εργασία σας↩

  8. Και εγώ στη θέση σας, το ίδιο θα έκανα↩

  9. Γιατί;↩

  10. Γιατί όμως να μην αλλάξουμε την ανάλυση των περιβαλλοντικών δεδομένων;↩

  11. Πώς τα αφαιρούμε από την μνήμη της R;↩

  12. Αυτό είναι το λεγόμενο hard boundary - άλλοι ερευνητές εφαρμόζουν το όριο \(> 0.8\) ή ακόμα και \(> 0.9\), αλλά είναι επιστημονικά ασφαλέστερο να εφαρμόζει κανείς το αυστηρότερο όριο↩

  13. Εδώ εφαρμόζουμε το αυστηρότερο όριο - μπορείτε να το αλλάξετε, εάν επιθυμείτε↩

  14. Σε απόσταση \(<1000m\) όταν εργαζόμαστε με χάρτες υποβάθρου κλίμακας 1 x 1 km2 (ή ανάλυσης 30 arc-seconds)↩

  15. Θα σας βοηθήσει το vignette του εν λόγου πακέτου↩

  16. Προτρέξτε στους Guisan et al. (2017) για περισσότερες λεπτομέρειες - εάν χρειάζεστε περισσότερη βιβλιογραφία, απευθυνθείτε σε εμένα (εάν φυσικά το επιθυμείτε)↩

  17. Όπως θα διαπιστώσετε παρακάτω, θα εφαρμόσουμε μια τεχνική η οποία ονομάζεται ensemble modelling. Σύμφωνα με αυτή, είναι προτιμότερο και επιστημονικά ορθότερο να μην χρησιμοποιεί κανείς ένα μόνο μοντέλο (π.χ., GLM ή τo MaxEnt)↩

  18. Προσέξτε ότι πρόκειται για αρχείο xlsx και όχι csv - αυτό σημαίνει ότι άνοιξα στο excel το αρχείο csv, χρησιμοποίησα την εντολή data to columns και το έσωσα ως αρχείο xlsx↩

  19. Δεν έχει σημασία ποιό θα διαλέξουμε - έχουν όλα την ίδια πιθανότητα↩

  20. Θα μπορούσε κάποιος να αντιτείνει ότι μπορούμε να χρησιμοποιήσουμε το πακέτο ecospat όταν \[\frac{\sum{Ανεξάρτητες,\, μη\, συγγραμικές\, μεταβλητές}}{\sum{Θέσεις\, εμφάνισης}} \geq 0.1\]↩

  21. Στο βιβλίο των Guisan et al. (2017) αναφέρεται λανθασμένα ότι ο αριθμός των pseudoabsences πρέπει να είναι τουλάχιστον 10000 - ο αριθμός αυτός έχει να κάνει με το MaxEnt, το οποίο ανήκει στην ομάδα των P/O μοντέλων↩

  22. Καλό είναι να είναι τουλάχιστον 10 τα set αυτά↩

  23. Θυμηθείτε να με ρωτήσετε τι σημαίνουν οι υπόλοιπες μεταβλητές↩

  24. Διαφορετικά δεν θα γίνει δεκτή και δεν θα βαθμολογηθείτε↩

  25. Πάντα σύμφωνα με το είδος σας↩

  26. Παραθέστε ΟΛΑ τα στοιχεία από τον σχετικό πίνακα↩

Κώστας Κουγιουμουτζής

Spring semester 2019-2020

 

Delivered to you by Dr. Kostas Kougioumoutzis

kkougiou@aua.gr